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CHAPTER  1 
INTRODUCTION 

Stress  wave  propagation  has  long  been  a  useful  tool  for  studying 
the  mechanical  properties  of  various  media,  especially  in  the  field 
of  material  science,  geophysics,  and  engineering.  Most  applications 
in  geotechnical  engineering  are  oriented^to  the  exploration  of  soil 
deposits  and  the  determination  of  soil  moduli  by  field  and  laboratory 
measurements  of  stress  wave  velocities.  These  soil  moduli  are  used 
in  dynamic  soil  structure  interaction  analyses  ranging  from  small- 
strain  problems,  as  associated  with  the  vibrations  of  machine  foundations, 
to  larger  strain  problems,  such  as  earthquake  or  blast  loading.  To 
simulate  the  situation  in  the  field,  knowledge  of  the  dynamic 
properties  of  soils  under  different  states  of  stress  -  isotropic, 
biaxial,  and  triaxial  -  is  important.  Unfortunately,  most  of  the 
research  results  available  have  emphasized  isotropic  stress  effects. 

A  research  program  was,  therefore,  undertaken  to  study  the  effects  of 
the  general  state  of  stress  on  shear  and  compression  wave  velocities  of 
soil  at  low  levels  of  strain.  The  studies  were  conducted  in  a  large 
scale  triaxial  testing  device  which  is  a  cube  constructed  of  reinforced 
stee^l,  designed  to  hold  a  cubic  soil  sample  measuring  7  ft  (2.1  m) 
on  each  side.  Rubber  loading  membranes  are  oriented  along  the  three 
principal  planes,  between  the  soil  and  three  mutually  perpendicular 
walls  of  the  cube.  The  membranes  are  hydraulically  pressurized  to 
attain  an  isotropic,  biaxial,  or  triaxial  state  of  stress.  Compression 
or  shear  waves  are  generated  along  the  principal  axes  through  one  of 


three  excitation  ports  located  on  the  three  walls  of  the  cube  which  are 
not  covered  with  membranes.  Accelerometer  packages  are  embedded 
in  the  soil  along  each  principal  axis  of  the  cube  to  monitor  the 
waves.  Soil  samples  used  to  date  have  been  dry  sand.  Details  and 
results  of  the  tests  have  been  presented  in  Kopperman,  et  al  (1982) 
and  Knox,  et  al  (1982) . 

While  most  of  the  effort  during  the  first  stage  of  the  research 
project  has  been  devoted  to  the  design  and  construction  of  the 
triaxial  cube,  some  work  has  also  been  performed  implementing  computer 
capabilities  which  could  be  used  in  the  future  for  interpretation  of 
the  experimental  data. 

The  first  problem  at  hand  was  to  determine  the  effect  of  th  e 
much  stiffer  (almost  rigid)  accelerometer  packages  embedded  in  the  sand 
on  the  characteristics  of  P-wave  and  S-wave  propagating  through  the 
soil  (scattering  or  filtering  effect).  Most  of  the  effort  in  this  work 
concentrated  on  shear  waves  because  of  the  greater  difficulty 
involved  in  these  measurements.  Two  approaches  were  considered  for 
the  solution  of  this  problem:  a  formulation  based  on  the  use  of  the 
boundary  integral  equation  (or  boundary  element)  method  and  a  more 
classical  discretization  using  finite  elements. 

The  boundary  integral  equation  method  is  particularly  appropriate 
to  study  the  filtering  effect  of  one  or  more  rigid  inclusions  on 
waves  propagating  through  a  full  elastic  space,  since  then  only  the 
boundaries  of  the  inclusions  must  be  discretized.  For  the  particular 
case  at  hand,  with  three  rigid  inclusions  along  any  principal  axis 
(and  considering  only  one  principal  axis),  formulation  of  the 


problem  in  the  frequency  domain  requires  the  solution  of  a  system  of 
54  equations  with  54  unknowns  for  each  frequency.  The  solution  for 
a  transient  excitation  can  then  be  obtained  using  Fourier  transforms. 

For  the  actual  case  of  the  cube,  by  opposition  to  the  full  space, 
it  would  also  be  necessary  to  discretize  its  outer  surface.  A 
discretization  appropriate  for  wavelengths  of  the  order  of  1  ft 
would  lead  to  a  system  of  56,502  equations  with  56,502  unknowns  for 
each  frequency.  By  taking  advantage  of  symmetry  and  antisymmetry 
conditions,  the  number  of  degrees  of  freedom  can  be  reduced  to  23,574, 
but  the  two  matrices  associated  with  tractions  and  displacements  are 
fully  populated,  which  makes  the  solution  prohibitively  expensive. 

If  three-dimensional  brick  elements  (with  linear  displacement 
expansion)  are  used,  an  appropriate  discretization  of  the  cube  into 
finite  elements  for  wavelengths  of  the  order  of  1  ft  or  larger 
requires  a  mesh  size  of  1.5  inches  or  56  elements  per  side  (1.5  inches 
is  also  the  size  of  the  rigid  inclusions).  This  leads  to  a  total  of 
175,616  elements  or  185,193  joints,  with  three  degrees  of  freedom  per 
joint  and  a  bandwidth  of  9582.  Solution  of  this  problem  would  again 
be  impractical  by  any  method  which  requires  assembling  the  total 
stiffness  matrix  or  solving  the  coupled  system  of  equations.  It  is 
possible,  on  the  other  hand,  to  solve  it  using  an  explicit  integration 
scheme,  such  as  the  central  difference  formula,  and  marching  out  the 
solution  in  the  time  domain.  While  the  time  step  of  integration  must 
be  kept  sufficiently  small  to  guarantee  stability  and  the  solution  is 
time  consuming,  the  advantage  of  this  approach  is  that  it  easily 
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allows  incorporation  of  nonhomo geneous  soil  properties  and  evantually 
nonlinear  soil  behavior. 

Because  of  these  considerations,  it  was  decided  that  the  finite 
element  formulation  had  definite  advantages  not  only  for  the  solution 
of  the  problem  at  hand  but  more  importantly  for  future  applications. 

By  taking  advantage  of  symmetry  and  antisymmetry,  one  fourth  of  the 
cube  was  discretized  using  28  x  28  x  55  elements  with  dimensions  of 
1.5  x  1.5  x  1.5  inches.  This  makes  the  side  of  the  cube  equal  to 
6  ft  10.5  inches  instead  of  the  7  ft  assumed,  but  this  small  difference 
allows  each  one  of  the  inclusions  to  be  represented  by  a  single  rigid 
element.  Even  with  these  simpl if actions ,  the  memory  capacity  of  the 
computer  available  at  the  University  of  Texas  (CDC  Dual  Cyber  170/750) 
did  not  allow  all  degrees  of  freedom  (141,288)  to  be  stored  in  core. 

A  computer  program  was  implemented  keeping  in  core  at  any  time  the 
degrees  of  freedom  corresponding  to  two  planes  (2523  x  2  =  5046)  and 
writing  on  disk  those  corresponding  to  one  plane  when  the  corresponding 
slice  of  finite  elements  had  been  processed,  while  reading  in  those  of 
a  new  plane.  Preliminary  studies  conducted  to  determine  the  smallest 
natural  period  and  the  required  time  step  of  integration  indicated  that 
some  90  to  100  steps  were  necessary  to  reproduce  the  arrival  of  the 
first  shear  wave  to  the  central  box.  For  this  duration  of  response, 
the  solution  of  the  three-dimensional  problem  takes  about  three  hours 
in  the  computer  of  the  University  of  Texas,  50  percent  of  this  time 
being  associated  with  the  transfer  of  data  from  primary  to  secondary 


memory  and  back. 


Since  only  a  small  number  of  fully  three-dimensional  analyses 
were  possible  within  the  budget  of  the  project,  it  was  decided  to 
implement  also  a  two-dimensional  (2-D) ,  plane-strain  solution  which 
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could  be  used  for  more  extensive  parametric  studies.  Two  versions  of 
this  program  were  written:  one  taking  again  advantage  of  anti-symmetry 
and  considering  a  mesh  of  55  x  28  elements  (3248  degrees  of  freedom) , 
and  the  other  with  the  complete  square  (55  x  55  elements  or  5272 
degrees  of  freedom).  The  smaller  was  used  to  verify  the  adequacy  of 
the  boundary  conditions  imposed  on  the  larger.  The  two-dimensional 
finite  element  model  is  shown  in  Fig.  1.1. 

Preliminary  analyses  were  performed  with  the  following  conditions: 
all  boundaries  were  fixed  except  for  the  small  window  area  on  the  left 
face  where  a  sinusoidal  acceleration  was  specified  with  a  frequency 
of  650  rad/sec  corresponding  to  a  wavelength  of  about  1  ft.  The  rigid 
inclusions  had  the  same  unit  mass  as  the  sand  sample.  The  results 
of  the  studies,  which  were  summarized  in  Stokoe,  et  al  (1982),  indicated 
that  for  the  range  of  wavelengths  of  interest,  while  there  is  a  small 
filtering  effect,  none  of  the  basic  characteristics  of  the  wave  (arrival 
time  of  the  first  shear  wave  and  its  shape)  would  be  significantly 
affected  by  the  presence  of  rigid  inclusions  of  the  size  of  the  accelero¬ 
meter  boxes.  Because  both  compression  waves  and  shear  waves  are 
generated,  the  first  shear  wave  arrival  time  is  difficult  to  estimate 
accurately.  It  is  easier  to  estimate  the  shear  wave  velocity  from 
the  time  interval  between  arrivals  of  the  waves  at  the  various  inclusions 
rather  than  from  the  time  of  first  arrival  at  one  location. 


Excitation 
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In  all  cases  studied  the  results,  in  the  form  of  accelerations 
versus  time,  showed  some  small  motion  in  the  direction  associated  with 
S-wave  motion  occurring  substantially  before  the  theoretical  time  of 
arrival  of  the  S-wave  (approximately  at  the  P-wave  arrival  time) . 

This  is  due  primarily  to  the  fact  that  the  excitation  is  applied  over 
a  small  area,  simulating  more  closely  a  point  source  and  a  spherical  wave 
front  (cylindrical  for  the  two-dimensional  case)  rather  than  a  pure  plane 
shear  wave  front.  On  the  other  hand  the  discrete  finite  element  model 
and  the  numerical  integration  scheme  will  clearly  introduce  some  errors, 
hopefully  small,  but  difficult  to  quantify.  The  boundary  conditions 
assumed  on  the  faces  of  the  cube  may  also  have  some  effect  on  the  type 
of  waves  that  reach  the  accelerometer  box. 

It  was  decided  therefore  to  investigate  further  these  effects  by 
considering  different  boundary  conditions,  different  types  of  excitation 
and  different  wave  lengths  (or  frequencies).  In  addition  to  the 
acceleration  time  histories  considered  in  the  initial  studies,  velocities 
and  displacements  were  plotted  in  order  to  assess  which  one  of  these 
variables  provided  an  easier  means  to  identify  initial  wave  arrival.  The 
results  of  these  additional  studies,  conducted  with  the  two-dimensional, 
plane  strain  model,  for  reasons  of  economy,  are  presented  in  this  report. 

Chapter  2  presents  the  results  relative  to  the  effect  of  the 
boundary  conditions  applied  on  the  faces  of  the  cube.  The  effect  of 
the  type  of  excitation  is  studied  in  Chapter  3,  and  the  effect  of  the 
frequency  excitation,  or  wavelength,  is  reported  in  Chapter  4.  Chapter  5 
contains  the  conclusions  of  the  study  and  recommendations  for  future  work. 


CHAPTER  2 


EFFECT  OF  BOUNDARY  CONDITIONS 


2.1  INTRODUCTION 

Two  types  of  boundary  conditions  are  considered  in  this  chapter: 
(1)  fixed  boundaries,  where  all  boundaries  are  assumed  to  be  fixed 
except  for  the  two  central  nodes  (for  2-D  case)  on  the  left  face  of 
the  cube,  where  a  sinusoidal  force  (instead  of  an  acceleration  as  in 
the  previous  report)  is  specified;  and  (2)  free  boundaries,  where  all 
boundaries  are  free  to  move  and  free  of  tractions  (except,  of  course, 
for  the  point  of  excitation).  Unfortunately,  for  this  second  case, 
the  four  boundaries  can  not  all  be  free  because  the  whole  system 
could  have  a  rigid  body  motion  causing  numerical  problems.  It  was 
decided  therefore  to  fix  the  boundary  on  the  right  side  (far  side  in 
the  path  of  wave  propagation)  because  most  of  the  waves  of  interest 
would  arrive  earlier  than  the  reflected  waves  from  this  boundary. 

The  frequency  of  excitation  used  for  these  studies  is  650  rad/sec 
corresponding  to  a  wavelength  of  approximately  1  ft.  The  rigid 
inclusions  and  soils  have  the  same  unit  mass  as  before. 

2.2  PRELIMINARY  CONSIDERATIONS 

All  results  in  the  last  report  were  presented  as  time  histories 
of  the  acceleration  at  three  selected  target  points.  These  target 
points  are  at  the  lower  left  corner  of  the  first,  second,  and  third 
inclusion  positions  which  are  referred  to  in  the  text  as  location 
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IX  *  14,  IX  =  28,  and  IX  =  42,  respectively  (see  Fig.  1.1).  In  this 
report,  some  results  at  the  point  of  excitation  will  also  be  shown  so 
that  changes  in  the  shape  of  the  waves  after  travelling  through  the 
medium  can  be  observed.  The  theoretical  arrival  times  for  plane  shear 
and  compressional  waves  to  each  target  point  are  indicated  in  the 
figures  by  vertical  lines  marked  with  the  letters  'S'  and  'P',  respectively. 

In  Fig.  2.1  the  displacement,  velocity  and  acceleration  for  a 
typical  case  are  plotted  on  the  same  time  scale.  It  can  be  seen, 
for  this  specific  case,  that  beyond  the  time  0.04  sec,  the  results 
behave  like  classical  harmonic  motion,  the  acceleration,  velocity 
and  displacement  being  nearly  sinusoidal  curves  with  a  constant  phase 
difference  between  each.  However,  the  earlier  motion  (before 
0.04  sec  in  this  case)  is  rather  complex  and  irregular,  starting  at  the 
theoretical  time  of  arrival  of  the  P-wave.  This  part  of  the  motion  is 
isolated  for  clarity  for  each  one  of  the  variables  (displacement, 
velocity  and  acceleration)  in  Fig.  2.2.  In  the  case  of  displacement, 
the  time  of  arrival  of  the  S-wave  is  only  slightly  ahead  of  the  zero 
crossing  point,  when  the  positive  motion  starts.  Points  "a",  corresponding 
to  the  zero  crossing,  or  Mb",  corresponding  to  the  time  where  the 
curve  starts  turning  upwards  are  often  used  in  experimental  work  to 
measure  the  arrival  time  of  the  S-wave.  If  the  same  approach  were 
followed  with  the  numerical  results,  the  estimated  shear  wave  velocity 
would  be  within  five  percent  of  the  exact  value. 

When  dealing  with  the  velocity  trace  on  the  other  hand,  if  point 
"d",  corresponding  to  the  zero  crossing,  were  chosen  the  estimated 
shear  wave  velocity  would  be  too  high,  with  an  error  of  10  to  15  percent. 
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Fig.  2.1  -  Displacement,  Velocity,  and  Acceleration  at  IX=28 
Free  Boundaries.  No  Rigid  Inclusions. 
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Fig.  2.2  -  Details  of  the  First  Arrival  of  Shear  Wave  from 
Displacement,  Velocity,  and  Acceleration  Curves 


The  theoretical  arrival  time  tends  to  be  about  one  third  up  the 
ascending  branch,  coinciding  approximately  with  the  inflection  point 
in  this  branch.  This  seems  also  to  be  consistent  with  experimental 
procedures  used  to  pick  the  initial  arrival  time. 

In  the  acceleration  trace  the  time  of  arrival  of  the  S-wave  is 
higher  up  on  the  ascending  branch  and  its  best  estimate  would  seem  to 
be  the  time  at  which  the  first  significant  positive  peak  occurs 
(point  "e") .  Use  of  this  point  would  yield  a  value  of  the  wave 
velocity  slightly  too  small,  with  an  error  of  about  5  percent. 

For  the  purposes  of  comparison,  Fig.  2.3  shows  the  displacement, 
velocity  and  acceleration  histories  at  the  point  of  application  of 
the  sinusoidal  force.  It  can  be  seen  that  the  displacement  starts 
with  a  slight  delay  due  to  the  numerical  solution  (no  correction 
was  applied  in  the  first  time  increment;  thus,  since  the  acceleration 
and  velocity  are  zero  at  t  =  0,  the  displacement  is  still  zero  at 
t  =  At).  The  time  at  which  the  displacement  curve  starts  coincides 
approximately  with  the  first  peak  in  the  acceleration  curve  and  the 
inflection  point  in  the  ascending  branch  of  the  velocity  trace. 

These  results  tend  to  indicate  that  while  there  are  some  errors 
associated  with  the  finite  element  discretization  and  the  numerical 
integration  scheme,  the  shear  wave  velocity  can  be  estimated  from 
the  numerical  results  with  reasonable  accuracy  (within  5  percent) 
if  the  appropriate  points  are  selected  in  the  time  histories. 

These  points  will  be  different  for  the  acceleration,  velocity  and 
displacement  traces, and  it  would  appear  that  the  displacement  curve 
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Fig.  2.3  -  Displacement,  Velocity,  and  Acceleration  at  the 
Point  of  Excitation.  Free  Boundaries.  No 
Rigid  Inclusions. 


14 


is  the  one  that  provides  the  clearest  reference  point.  This 
conclusion  seems  logical  if  one  takes  into  account  that  the  numerical 
integration  procedure  used  in  the  analyses  computes  directly  the 
values  of  the  displacements  versus  time,  and  that  velocities  and 
accelerations  must  then  be  obtained  from  the  displacements.  Since 
the  velocity  and  accelerations  are  the  first  and  second  derivatives 
of  the  displacements,  these  time  histories  are  less  smooth. 

2.3  EFFECT  OF  BOUNDARY  CONDITIONS  ON  S-WAVES 

Figures  2.4  to  2.6  show  the  time  histories  of  displacements, 
velocities  and  accelerations  at  one  of  the  points  of  application  of 
the  load  for  the  cases  of  fixed  and  free  boundaries  and  no  rigid 
inclusions  within  the  soil  mass.  In  both  cases  the  displacement  curve 
starts  slightly  after  the  zero  time.  As  pointed  out  earlier,  this  is 
due  to  the  starting  procedure  in  the  numerical  algorithm  which  forces 
a  zero  displacement  at  time  At.  This  could  be  improved  by  applying 
a  correction  in  the  first  step  of  integration  (typically  using  a 
different  formula  for  this  step).  All  three  figures  show  that  the 
amplitudes  of  motion  (displacement,  velocity  or  acceleration)  are 
similar  for  both  types  of  boundaries.  The  results  for  the  free 
boundaries  are,  however,  smoother,  which  can  be  seen  in  the  velocity 
and  particularly  in  the  acceleration  traces.  They  also  tend  to  lag 
behind  the  results  of  the  fixed  case. 

Figures  2.7  to  2.9  and  2.10  to  2.12  show  the  corresponding  results 
at  the  target  points  IX  =  14  and  IX  =  28  (left  lower  corners  of  the 
positions  of  the  first  and  second  inclusions),  respectively.  The  same 
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Fig.  2.4  -  Vertical  Displacements  at  the  Point  of  Excitation 
Fixed  and  Free  Boundaries.  No  Rigid  Inclusions. 
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Fig.  2.6  -  Vertical  Accelerations  at  the  Point  of  Excitation 
Fixed  and  Free  Boundaries.  No  Rigid  Inclusions. 
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Fig.  2.8  -  Vertical  Velocities  at  IX=14.  Fixed  and  Free 
Boundaries.  No  Rigid  Inclusions. 
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observations  can  be  made,  but  in  these  cases,  in  addition,  the  motion 
for  the  fixed  boundaries  has  clearly  a  smaller  amplitude  than  that  of 
the  free  model. 

The  corresponding  results  with  three  rigid  inclusions  are  presented 
in  Figures  2.13  to  2.15  (motion  at  excitation  point),  2.16  to  2.18 
(motion  at  IX  =  14)  and  2.19  to  2.21  (motion  at  IX  =  28).  The  presence 
of  these  rigid  inclusions  does  not  produce  any  significant  difference 
in  the  shape  of  the  waves  or  in  the  effect  of  the  boundary  conditions, 
for  the  size  of  the  inclusions  and  the  wavelength  considered. 

Since  the  results  for  the  free  boundaries  lag  in  all  cases  behind 
those  corresponding  to  fixed  boundaries,  the  shear  wave  velocity 
estimates  based  on  the  first  arrival  time  would  be  smaller  for  the 
former  case  than  for  the  latter.  On  the  other  hand,  when  interval 
arrival  times  between  targets  are  used,  by  picking  the  corresponding 
points  on  the  curves  for  each  target  position,  the  calculated  shear  wave 
velocities  are  about  the  same,  irrespective  of  the  boundary  conditions. 
Figures  2.22  through  2.27  show  interval  arrival  times  when  there  are 
no  inclusions  and  Figs.  2.28  through  2.33  show  the  corresponding  results 
with  three  rigid  inclusions.  When  the  interval  times  are  measured 
between  the  first  S-wave  peaks,  calculated  shear  wave  velocities  have 
average  values  of  101.2  ft/sec  for  free  boundaries  and  no  inclusions, 
102.4  ft/sec  for  free  boundaries  and  three  inclusions,  and  97  and 
97.5  ft/sec  for  fixed  boundaries  (no  inclusions  and  three  inclusions, 
respectively).  The  values  for  free  boundaries  tend  to  be,  therefore, 

slightly  higher  than  those  for  fixed  boundaries,  by  opposition  to  the 

estimates  based  in  initial  arrival  of  the  S-wave.  Measuring 

interval  travel  times  from  the  second  peaks,  a  less  desirable 
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Fig.  2.13  -  Vertical  Displacements  at  the  Point  of  Excitation. 

Fixed  and  Free  Boundaries.  Three  Rigid  Inclusions 


Fig.  2.15  -  Vertical  Accelerations  at  the  Point  of  Excitation. 

Fixed  and  Free  Boundaries.  Three  Rigid  Inclusions 
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Fig.  2.16  -  Vertical  Displacements  at  IX=14.  Fixed  and 
Free  Boundaries.  Three  Rigid  Inclusions. 
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Fig.  2.19  -  Vertical  Displacements  at  IX=28.  Fixed  and 
Free  Boundaries.  Three  Rigid  Inclusions. 
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solution,  since  it  involves  a  longer  time  of  integration,  the  corres¬ 
ponding  values  are  96.7,  97.7,  96  and  96.5  ft/sec,  respectively. 

These  values  are  consistently  smaller  and  less  accurate  (the  actual 
shear  wave  velocity  of  the  soil  used  in  the  analyses  was  100  feet  per 
second) .  It  can  also  be  noticed  that  the  presence  of  the  inclusions 
increases  the  value  of  the  calculated  shear  wave  velocity  by  about 
one  percent  at  most. 

Figures  2.34  to  2.36  illustrate  better  the  effect  of  the  rigid 
inclusions  for  the  case  of  free  boundaries  and  Figs.  2.37  to  2.39 
for  fixed  boundaries.  It  can  be  seen  that  the  shapes  are  practically 
the  same  for  no  inclusions,  one  inclusi  n  or  three  inclusions.  Changes  in 
the  arrival  times  or  position  of  the  peaks  are  very  small  leading  to 
differences  in  the  estimated  wave  velocites  of  the  order  of  one  percent 
The  main  effect  of  the  inclusions  is  a  smoothing  of  the  curves,  particu¬ 
larly  noticeable  in  the  acceleration  traces  of  the  fixed  boundary 
case,  and  a  small  reduction  in  amplitude  (filtering).  There  is  a 
reduction  with  a  placement  of  a  single  inclusion  at  the  center  and 
an  additional  one  with  the  placement  of  the  other  two  inclusions.  The 
filtering  effect  seems  to  be  of  the  same  nature  and  magnitude  for 
both  sets  of  boundary  conditions. 

2.4  EFFECT  OF  BOUNDARY  CONDITIONS  ON  P-WAVES 

Most  of  the  studies  conducted  assumed  a  vertical,  shear  type 
excitation,  since  it  was  felt  that  determination  of  the  shear  wave 
velocity  might  present  more  problems  than  that  of  the  P-wave.  A 
number  of  cases  were  considered,  however,  with  a  horizontal,  compressional 
type  force  which  would  generate  primarily  P-waves. 
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Fig.  2.38  -  Free  Field  and  Inclusion  Velocities.  Fixed  Boundaries 
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Fig.  2.39  -  Free  Field  and  Inclusion 
Boundaries. 


Figures  2.40  to  2.42  show  the  time  histories  of  the  horizontal 
displacement,  velocity  and  acceleration  at  target  point  1X28  (left 
corner  of  the  middle  inclusion),  when  there  are  no  inclusions,  for 
free  and  fixed  boundaries.  The  motion  starts  now  only  slightly  ahead 
of  the  theoretical  arrival  time,  with  a  well  defined  beginning. 

All  three  (or  six)  curves  start  with  a  small  positive  curvature 
which  can  be  attributed  to  the  integration  scheme.  If  the  displacement 
curves  were  rounded  at  the  start  to  make  them  approach  more  closely 
a  half  sinusoid,  the  theoretical  arrival  time  of  the  P-wave  would  be 
predicted  very  well.  For  the  velocity  curves,  on  the  other  hand, 
the  arrival  time  would  correspond  again  more  closely  to  the  inflection 
point  on  the  ascending  branch,  while  for  the  acceleration  it  is  just 
slightly  ahead  of  the  first  peak.  As  in  the  case  of  S-waves,  the 
first  time  of  arrival  would  be  easier  to  estimate  from  the  displacement 
traces,  which  are  smoother  curves. 

The  motions  for  the  model  with  free  boundaries  lag  again  behind 
those  for  fixed  boundaries,  except  at  the  start  where  they  are  almost 
identical.  On  the  other  hand  the  amplitudes  of  motion  are  larger 
for  fixed  boundaries,  contrary  to  what  was  observed  for  a  shear  type 
excitation. 

The  effect  of  the  rigid  inclusions  is  illustrated  in  Figs.  2.43 
to  2.45,  where  the  results  for  the  cases  of  no  inclusions,  one  inclusion 
and  three  inclusions  are  directly  compared  for  the  solution  with  free 
boundaries.  The  main  effect  is  again  a  filtering  of  higher  frequency 
components,  particularly  in  the  acceleration  trace  (Fig.  2.45)  and 
a  reduction  in  the  amplitudes  of  motion. 


Fig.  2.40  -  Horizontal  Displacements  at  IX_28.  Fixed  and  Free 
Boundaries.  No  Rigid  Inclusions.  P-wave. 
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Fig.  2.42  -  Horizontal  Accelerations  at  IX*28.  Fixed  and 
Boundaries.  No  Rigid  Inclusions.  P-wave. 
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CHAPTER  3 

EFFECT  OF  TYPE  OF  EXCITATION 


3.1  INTRODUCTION 

As  mentioned  in  Chapter  1,  most  studies  were  conducted  attempting 
to  simulate  the  conditions  existing  in  the  cube  during  the  experimental 
work.  As  a  result  the  loads  were  applied  over  a  small  area,  simulating 
more  closely  a  point  source  and  a  spherical  wave  front  (cylindrical 
for  the  two-dimensional  model)  than  a  plane  wave.  Under  these 
conditions  vertical  motions  due  to  a  vertical,  shear  type  excitation, 
were  starting  at  the  time  of  arrival  of  the  P-waves.  To  determine  to 
what  extent  this  was  due  to  the  type  of  excitation,  and  reproduced 
therefore  the  physical  reality,  or  to  noise  introduced  by  the 
discretization  and  the  integration  scheme,  additional  studies 
were  conducted  changing  the  type  of  excitation.  Three  different 
situations  were  considered:  one  in  which  the  point  excitation  was 
replaced  by  equal  forces  at  all  points  on  the  left  boundary  simulating 
a  uniformly  distributed  harmonic  force;  a  second  one  with  still  a 
point  load  but  preventing  the  horizontal  displacement  of  all  nodes; 
and  a  third  one  with  the  distributed  load  and  the  same  constraints  on 
the  nodal  displacements.  The  results  for  each  one  of  these  models 
are  discussed  in  the  following  sections.  In  all  cases  the  free  boundary 
conditions  were  applied, and  it  was  assumed  that  there  were  no  inclusions. 

3.2  UNIFORMLY  DISTRIBUTED  EXCITATION 

The  model  used  in  this  case  was  the  same  discussed  in  Chapter  2 
for  free  boundaries,  the  only  difference  being  in  the  applied  loads. 


The  nodal  forces  were  scaled  down  so  as  to  have  the  same  magnitude  of 
total  load. 

Figures  3.1  and  3.2  compare  the  time  histories  of  accelerations 
at  the  target  points  1X14  and  1X28  for  the  point  load  and  the  distributed 
force.  It  can  be  seen  that  with  the  distributed  excitation  the  small 
motion  that  started  at  the  time  of  arrival  of  the  P-wave  no  longer 
exists.  The  motion  starts  still  a  little  early,  with  a  positive 
curvature,  but  if  the  curve  were  to  be  rounded  to  approximate  a 
half  sinusoid,  the  starting  point  would  almost  coincide  with  the 
theoretical  time  of  arrival  of  the  S-wave  (in  spite  of  the  fact  that 
these  are  acceleration  traces) .  The  motion  is  nearly  a  perfect 
sinusoid  with  constant  amplitude,  particularly  for  the  second  target 
point.  The  amplitude  of  motion  is  smaller  than  for  the  case  of  a 
concentrated  load,  as  could  be  expected,  and  the  peaks  are  slightly  out 
of  phase,  occurring  later  for  the  distributed  load. 

3.3  POINT  EXCITATION  WITH  NO  HORIZONTAL  DISPLACEMENTS 

In  this  second  case  the  point  excitation  is  again  used  but  all 
nodes  are  forced  to  move  in  vertical  planes  (vertical  lines  actually) 
to  simulate  the  kind  of  motion  that  would  take  place  under  a  plane 
shear  wave  travelling  in  the  horizontal  direction. 

Figures  3.3  and  3.4  show  again  the  vertical  accelerations  at  the 
target  points  1X14  and  1X28,  for  this  case  and  the  standard  solution 
of  Chapter  2.  It  can  be  seen  that  when  horizontal  displacements  are 
prevented  the  motion  no  longer  starts  at  the  time  of  arrival  of  the 
P-waves.  It  starts,  however,  substantially  before  the  theoretical 
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S-wave  arrival  time.  For  the  first  few  cycles  the  motion  is  nearly  in 
phase  with  the  results  without  displacement  constraints.  As  a  conse¬ 
quence  the  theoretical  time  of  arrival  of  the  S-wave  is  still  relatively 
close  to  the  first  peak.  For  later  cycles  the  constrained  solution  is 
slightly  ahead  of  the  standard  one,  but  their  amplitudes  tend  to  be 
more  similar. 

Preventing  horizontal  displacements  but  applying  still  a  concen¬ 
trated  force  does  not  eliminate  enirely  volumetric  deformations.  The 
motions  are  not  as  smooth  and  regular  as  for  the  case  of  a  distributed 
force. 

3.4  UNIFORM  LOAD  WITH  NO  HORIZONTAL  DISPLACEMENTS 

The  changes  introduced  in  the  model  in  the  previous  two  sections 
were  combined,  using  a  uniformly  distributed  excitation  on  the  left 
face  (instead  of  a  point  force)  and  preventing  at  the  same  time  the 
hroizontal  displacements  of  all  the  nodes.  The  reuslts  are  shown  in 
Figs.  3.5  and  3.6  and  are  almost  identical  to  those  obtained  in 
Section  3.2  (uniform  load  but  no  displacement  constraints).  The  motion 
resembles  again  a  sinusoidal  curve  with  uniform  amplitude  and  starts 
only  slightly  ahead  of  the  theoretical  time.  Figure  3.7  compares  the 
results  at  1X28  with  those  of  Section  3.2  indicating  that  the 
enforcement  of  only  vertical  displacements  has  very  little  effect  when 
the  applied  loads  are  uniformly  distributed  over  the  left  face  (it 
should  have  no  effect  at  all  for  a  pure  plane  shear  wave  travelling 
in  the  horizontal  direction) . 


Vertical  Acceleration,  ft/sec 
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Fig.  3.5  -  Vertical  Accelerations  at  IX*14.  Point  Excitation  and 
Uniformly  Distributed  Excitation  with  Restraints  on 
Horizc  tal  Movement. 
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CHAPTER  4 


EFFECT  OF  WAVELENGTH 

4.1  INTRODUCTION 

All  studies  described  in  the  previous  chapters  and  the  first 
report  were  based  on  a  sinusoidal  excitation  with  a  frequency  of 
650  rad/sec,  corresponding  to  a  wavelength  of  about  one  foot  (1.03  ft 
actually) .  This  value  was  selected  because  the  wavelength  is  approximately 
eight  times  the  size  of  the  inclusion,  and  it  was  felt  it  would  be  the 
shortest  wavelength  used  in  the  experiments.  The  effect  of  the  inclu¬ 
sions  should  be  even  smaller  for  longer  wavelengths  (smaller  frequencies) 
but  the  effect  of  the  boundaries  (finite  dimensions  of  the  cube)  could 
be  more  pronounced.  The  size  of  the  inclusions  is  also  the  size  of 
each  element  in  the  finite  element  mesh.  With  eight  elements  per 
wavelength,  the  numerical  results  should  still  have  good  accuracy.  It 
is  commonly  accepted  that  in  dynamic  analyses  with  fi  ite  elements  the 
element  size  should  be  smaller  than  1/4  to  1/6  of  the  wavelength  to 
obtain  reasonable  answers.  This  implies  that  somewhat  higher  frequencies 
could  be  used  with  the  model  implemented,  but  a  frequency  of  1250  rad/sec 
would  be  a  practical  upper  limit. 

Additional  analyses  were  conducted  using  two  different  frequencies 
of  excitation:  a  frequency  of  325  rad/sec,  corresponding  to  a  wave¬ 
length  of  approximately  2  ft  (or  16  times  the  inclusion  size),  and  a 
frequency  of  975  rad/sec,  corresponding  to  a  wavelength  of  about  0.65  ft 
(about  five  times  the  size  of  the  inclussion).  In  this  chapter,  these 
frequencies  will  be  referred  to  as  'low'  and  'high'  frequencies,  respectively. 
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The  650  rad/sec  used  In  the  other  studies  will  be  considered  as  the 
'normal'  frequency.  For  these  new  analyses  only  the  case  with  free 
boundary  conditions  was  considered. 

4.2  AMPLITUDE  AND  SHAPE  OF  WAVES 

Figures  4.1  and  4.2  show  the  displacement,  velocity  and  acceleration 
records  at  the  point  of  excitation  for  the  low-frequency  and  high- 
frequency  excitations,  respectively.  These  results  are  directly 
comparable  to  those  in  Fig.  2.3  for  the  normal  frequency.  The  shapes 
and  the  general  characteristics  of  the  motions  are  similar,  and  the 
observations  made  in  Chapter  2  would  seem  to  apply  to  all  three 
frequencies.  The  acceleration  trace  for  the  low-frequency  excitation 
shows,  however,  higher  frequency  components,  which  are  particularly 
clear  in  the  first  cycles  and  appear  to  get  smoothed  out  in  the 
following  ones. 

The  displacement,  velocity  and  acceleration  at  the  central  target 
point  1X28  are  shown  in  Figs.  4.3  to  4.5  for  all  three  frequencies. 

The  amplitude  of  the  first  positive  displacement  peak  increases 
markedly  as  the  frequency  decreases.  This  is  still  true  for  the 
velocity  and  even  for  the  acceleration  although  the  increase  is  smaller 
in  each  case.  The  first  acceleration  peak  is  poorly  defined  for  the 
high-frequency  excitation.  For  the  following  peaks  the  amplitude 
of  the  displacement  seems  to  be  largest  for  the  normal  frequency  and 
smallest  for  the  high  frequency  but  there  are  not  enough  cycles  of 
the  low  frequency  response  for  a  meaningful  comparison  (the  displacement 
amplitude  at  low-frequency  could  in  fact  be  larger) .  The  normal 
frequency  has  clearly  the  largest  amplitude  of  motion  in  the  velocity 
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4.5  -  Vertical  Accelerations  at  IX=28  for  Different 
Frequencies  of  Excitation. 


trace,  and  the  low  frequency  the  smallest.  For  the  acceleration 
records  the  high  frequency  has  the  largest  amplitude,  followed  by  the 
normal  and  the  low  frequency.  This  is  proper  since  the  velocity  for 
a  steady  state  condition  would  be  equal  to  the  displacement  multiplied  by 
the  circular  frequency,  while  the  acceleration  would  be  the  displacement 
multiplied  by  the  square  of  the  frequency.  The  acceleration  trace  for 
the  low  frequency  at  1X28  shows  again  higher  frequency  components. 

The  theoretical  time  of  arrival  of  the  S-wave  coincides  almost 
exactly  with  the  zero  crossing  point  of  the  displacement  curve  for  the 
high  frequency,  is  only  slightly  ahead  of  this  point  for  the  normal 
frequency  and  is  further  ahead  for  the  low  frequency.  The  estimate 
of  the  shear  wave  velocity  based  on  the  zero  crossing  point  would  be 
about  five  percent  too  small  for  the  low  frequency.  This  shift  to  the 
right  of  the  curves  with  decreasing  frequency  is  also  apparent  in  the 
velocity  and  acceleration  curves.  The  finite  element  model  and  the 
integration  scheme  would  seem  to  result  in  a  dispersive  solution,  with 
the  measured  shear  wave  velocity  being  a  function  of  the  frequency  of 
excitation,  if  the  velocity  were  determined  from  the  estimate  of  the 
initial  arrival. 

To  explore  further  this  matter  Figs.  4.6  to  4.8  show  the  corre¬ 
sponding  results  at  the  third  target  point,  1X42.  The  sift  in  the  curves 
is  still  apparent,  but  the  theoretical  time  of  arrival  of  the  S-wave  is 
now  in  between  the  zero  crossing  points  for  the  high-  (and  normal)  and 
low-frequency  records.  The  shear  wave  velocity  estimated  from  the 
displacement  curves  (Fig.  4.6)  would  be  about  three  percent  too 
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Fig.  4.7  -  Vertical  Velocities  at  IX-42  for  Different 
Frequencies  of  Excitation. 
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8  -  Vertical  Accelerations  at  IX*A2  for  Different 
Frequencies  of  Excitation. 


large  for  the  high  and  normal  frequency  and  about  one  percent  (or  less) 
too  small  for  the  low  frequency.  The  same  would  be  basically  true 
for  the  velocity  curves  using  the  inflection  point  on  the  ascending 
branch  to  estimate  the  time  of  arrival.  With  the  acceleration  curve, 
using  the  first  peak  for  the  arrival  time,  determination  is  difficult 
for  the  high  frequency  excitation  because  the  peak  is  not  well  defined. 

For  the  normal  frequency  the  shear  wave  velocity  would  be  underestimated 
by  less  than  one  percent,  and  for  the  low  frequency  it  would  be  under¬ 
estimated  by  about  three  percent. 

Thus,  while  the  shift  in  the  curves  with  decreasing  frequency 
occurs  at  all  points  for  displacement,  velocity  and  acceleration, 
whether  it  results  in  an  underestimation  of  the  shear  wave  velocity 
at  low  frequencies  or  an  overestimation  at  high  frequencies  depends 
on  the  target  point  and  the  particular  curve.  For  the  range  of 
frequencies  considered,  which  is  similar  to  the  range  of  interest 
in  the  experimental  work,  the  errors  are,  however,  of  the  order  of 
five  percent  or  less. 

Figures  4.9  to  4.11  illustrate  the  detemrination  of  the  shear  wave 
velocity  using  interval  times  for  the  high-frequency  excitation  and 
Figs  4.12  to  4.14  for  the  low-frequency  case.  Using  the  first  peaks 
produces  again  better  results  than  using  the  second.  For  the  displace¬ 
ment  curves  and  the  high  frequency,  estimates  of  the  wave  velocity  are 
97.2  and  102.9  ft/sec  (first  to  second  target  points  and  second  to  third 
target  points,  respectively),  with  an  average  of  almost  exactly  100  ft/sec 
For  the  low  frequency  the  corresponding  values  are  104.8  and  101.2  ft/sec, 
with  an  average  of  103  ft/sec.  Notice  that  this  trend  is  opposite  to 
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the  one  discussed  above  for  first  arrival  times,  where  the  results  for 
the  low  frequency  excitation  would  yield  a  slightly  lower  velocity. 

By  using  the  velocity  time  histories,  the  estimates  of  wave  velocity  are 
100  and  102.9  fts/ec  (average  of  101.5)  for  the  high  frequency,  and  the 
same  for  the  low  frequency.  With  the  acceleration  curves,  no  estimate 
can  be  obtained  from  the  first  peaks  for  high  frequency  because  of 
their  poor  definition;  for  the  low  frequency  the  va  ues  are  100  and 
104.2  ft/sec,  with  an  average  of  102.1  ft/sec. 

4.3  EFFECT  OF  RIGID  INCLUSIONS 

Figures  4.15  to  4.17  compare  the  displacement,  velocity  and 
acceleration  time  histories  at  target  point  1X28  (left  corner  of  the 
position  of  the  cnetral  inclusion)  for  the  cases  of  no  inclusions, 
one  inclusion  and  three  inclusions  and  the  high-frequency  excitation. 
Figures  4.18  to  4.20  show  the  corresponding  results  for  the  low- 
frequency  case. 

The  presence  of  the  inclusions  produces  again  some  very  minor 
shifts  in  the  curves,  but  the  main  effect  is  the  filtering  of  the 
motion.  The  reduction  in  amplitude  is  significant  for  the  high- 
frequency  case  and  more  pronounced  than  for  the  normal-frequency  case 
(Figs.  2.34  to  2.39),  as  could  be  expected.  It  is  very  small  for 
low-frequency  excitation,  although  there  is  some  smoothing  of  the  high- 
frequency  components  in  the  acceleration  trace. 

The  effect  of  the  inclusions  on  the  initial-arrival  or  interval 
travel  times  and  on  the  estimates  of  the  shear  wave  velocity  is  very  small 
and  independent  of  the  frequency  of  the  excitation  (within  the  accuracy  of 
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Fig.  A. 16  -  Free  Field  and  Inclusion  Velocities 
Boundaries.  High  Frequency. 
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Fig.  4.17  -  Free  Field  and  Inclusion  Accelerations 
Boundaries.  High  Frequency. 


these  estimates  and  for  the  range  of  parameters  studied).  The  reduction 
in  the  amplitude  of  the  waves  increases,  as  expected,  as 
the  wavelength  decreases  and  approaches  the  size  of  the  inclusion, 
but  this  is  not  of  importance  for  the  studies  conducted  with  the 


CHAPTER  5 


CONCLUSIONS 


The  studies  presented  in  this  report  were  conducted  to  assess 
further  the  ability  of  a  finite  element  model  with  an  explicit  integration 
scheme  to  reproduce  the  propagation  of  waves  in  a  soil  cube  due  to  a 
harmonic  excitation  applied  over  a  small  area  on  one  of  the  faces.  Of 
particular  interest  was  the  study  of  the  effect  of  the  boundary  conditions 
applied  on  the  faces  of  the  cube,  the  type  of  excitation,  and  the  excitation 
frequency  or  wavelength,  on  the  shape  of  the  waves,  their  times  of  arrival 
and  on  estimates  of  wave  velocity  obtained  form  the  results.  Due 
to  the  cost  of  computation  the  analyses  were  performed  with  a  two- 
dimensional,  instead  of  a  fully  three-dimensional,  model.  It  is  believed, 
however,  that  the  basic  conclusions  derived  from  these  studies  are 
applicable  to  both  cases. 

The  results  presented  in  the  previous  chapters  indicate  that: 

1.  When  dealing  with  a  shear  type  excitation  applied  over  a  small 
area,  which  could  almost  be  considered  a  point  load,  the  corresponding 
motions  at  various  target  points  will  start  at  the  theoretical  time  of 
arrival  of  the  P-wave  and  not  at  the  arrival  time  of  the  S-wave.  This 
early  motion  is  more  significant  when  the  target  point  is  relatively 
close  to  the  point  of  excitation  and  tends  to  flatten  out  as  the 
distance  increases.  When  the  excitation  is  a  uniformly  distributed 
force,  simulating  more  closely  a  plane  wave  front,  this  early  motion 
does  not  occur,  and  motion  starts  only  slightly  ahead  of  the  theoretical 
time.  This  arrival  time  can  be  predicted  very  well  by  rounding  the  first 


half  of  the  wave. 
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2.  When  dealing  with  a  compressional  type  excitation  applied  over 
a  small  area,  the  above  problem  no  longer  exists,  but  the  motion  still 
starts  slightly  ahead  of  the  thoeretical  time. 

3.  For  a  point  type  excitation,  if  wave  propagation  velocity  is 
to  be  estimated  from  the  time  of  arrival  of  the  first  wave,  selection 
of  the  arrival  point  depends  on  the  type  of  trace  used,  displacement, 
velocity  or  acceleration.  The  arrival  point  should  be  selected  as 

the  zero  crossing  point  or  the  starting  point  of  the  first  wave  for  the 
displacement  curve,  the  inflection  point  on  the  ascending  branch  of  this 
wave  for  the  velocity  record,  and  a  point  close  to  the  peak  of  this 
wave  (and  perhaps  slightly  before  the  peak)  for  the  acceleration  trace. 
The  displacement  time-history  curve  is  the  smoothest  one  and  seems  to 
provide  the  easiest  reference  point.  The  velocity  curve  is  also  useful, 
but  the  acceleration  record  may  not  have  the  first  peak  well  defined, 
particularly  for  high  frequencies. 

4.  Determination  of  the  wave  propagation  velocity  from  the  interval 
travel  times  between  target  points  is  somewhat  easier  and  more  reliable 
than  from  the  first  arrival.  Taking  three  target  points  and  averaging 
the  values  resulting  from  the  travel  times  between  the  targets  seems 

to  provide  the  best  results. 

5.  Boundary  conditions  applied  on  the  faces  of  the  model  affect 
the  shape  and  patterns  of  the  waves  but  has  only  a  small  influence  on 
arrival  times.  The  motions  for  free  boundary  conditions  lag  slightly 
behind  those  obtained  with  fixed  boundaries.  This  implies  that  if  the 
wave  propagation  velocites  are  estimated  from  the  first  arrival  time 
they  would  be  somewhat  smaller  (about  five  percent)  for  the  free  boundary 


model.  On  the  other  hand  the  velocities  estimated  from  interarrival 
times  are  slightly  larger  with  free  boundary  conditions. 

6.  As  the  frequency  of  the  excitation  decreases,  within  the  range 
of  frequencies  studied,  the  arrival  of  the  first  wave  seems  to  be  delayed. 
Wave  velocities  determined  from  the  first  arrival  time  would  decrease 
slightly  with  decreasing  frequency  (increasing  wavelength).  As  in  item 

5  the  opposite  is  true  when  measuring  interval  travel  times. 

7.  The  presence  of  rigid  inclusions  of  the  dimensions  considered 
herein  (one  sixteenth  to  one  fifth  of  the  wavelength,  at  a  distance 
equal  to  fourteen  times  their  size)  has  only  a  very  minor  effect  on  the 
first-arrival  and  interval  times.  Estimated  shear  wave  velocites  with 
and  without  inclusions  change  by  less  than  one  percent.  The  inclusions 
smooth  slightly  the  motions  (particularly  the  acceleration  traces)  and 
reduce  the  amplitudes.  This  reduction  is  negligible  for  a  wavelength 
sixteen  times  the  inclusion  size,  small  for  half  this  wavelength  and 
more  significant  when  the  wavelength  is  only  five  times  the  size  of  the 
inclusion. 

8.  While  the  numerical  solution  has  some  errors  it  appears  that 

for  the  range  of  parameters  studied,  with  some  care,  shear  wave  (or  P-wave) 
velocities  could  be  predicted  from  the  results  of  the  analyses  within 
five  percent.  This  should  be  sufficient  accuracy  for  most  practical 
applications. 

It  was  intended  in  all  these  studies  to  simulate  the  wave  patterns 
that  would  be  generated  inside  the  cube  during  the  experimental  work. 

This  is  the  reason  for  the  finite  domain,  the  size  and  position  of  the 
Inclusions  and  the  range  of  frequencies  used  for  the  excitation.  The 


more  general  area  of  wave  propagation  in  a  nonhomogeneous  medium,  or  a 
medium  with  inclusions  of  different  properties  (as  might  be  an  alluvial 
deposit),  is  one  that  has  received  relatively  little  attention  and  which 
deserves  further  research.  By  considering  a  full  space,  and  thus  elimina 
ting  the  effect  of  the  boundary  conditions,  the  accuracy  of  the  finite 
element  solution  can  be  better  evaluated  by  comparing  the  results  to 
those  obtained  with  the  boundary  element  (or  boundary  integral  equation) 
method.  The  effects  of  size,  position  and  properties  of  the  inclusions 
for  a  more  extensive  range  of  frequencies  should  then  be  investigated. 
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